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^* ■ Abstract. The maximally compact representation of a regular orbit 

is in terms of its action-angle variables (J,0). Computing the map be- 
tween a trajectory's Cartesian coordinates and its action-angle variables 
is called torus construction. This article reviews various approaches to 
torus construction and their application to galactic dynamics. 
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1. Introduction 



In systems with a single degree of freedom, constancy of the energy allows the 
momentum variable p to be written in terms of the coordinate variable q as 

^O ' H[p, q) = E, and the dependence of both variables on time follows immediately 

^^ . from Hamilton's equations. In general systems with N > 2 degrees of freedom 

(DOF), such a solution is generally not possible unless the Hamilton- Jacobi 
equation is separable, in which case the separation constants are isolating inte- 
grals of the motion. An isolating integral is a conserved quantity that in some 
transformed coordinate system makes dH/dpi = f{qi), thus allowing the motion 
O I in gj to be reduced to quadratures. Each isolating integral restricts the dimen- 

4— > ' sionality of the phase space region accessible to an orbit by one; if there are N 

j^ I such integrals, the orbit moves in a phase space of dimension 2N — N = N, and 

the motion is regular. The A^-dimensional phase space region to which a regular 
orbit is confined is topologically a torus (Figure 1). Orbits in time-independent 

^ ' potentials may be either regular or chaotic, respecting a smaller number of in- 

tegrals - typically only the energy integral E. Although chaotic orbits are not, 
strictly speaking, confined to tori, numerical integrations suggest that many 
chaotic trajectories are effectively regular, remaining confined for long periods 
of time to regions of phase space much more restricted than the full energy 
hypersurface. 

The most compact representation of a regular orbit is in terms of the coor- 
dinates on the torus (Figure 1) - the action-angle variables (J, 6). The process of 
determining the map (x, v) -^ (J,^) is referred to as torus construction. There 
are a number of contexts in which it is useful to know the (J, ^). One example is 
the response of orbits to slow changes in the potential, which leave the actions 
(J) unchanged. Another is the behavior of weakly chaotic orbits, which may be 
approximated as regular orbits that slowly diffuse from one torus to another. 
A third example is galaxy modeling, where regular orbits are most efficiently 
represented and stored via the coordinates that define their tori. 

This article reviews techniques for mapping Cartesian coordinates into 
action-angle variables in non-integrable potentials. Two general approaches to 
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this problem have been developed. Trajectory- following algorithms are based 
on the quasi-periodicity of regular motion: Fourier decomposition of the tra- 
jectory yields the fundamental frequencies on the torus as well as the spectral 
amplitudes, which allow immediate construction of the map 6 ^> x. Iterative 
approaches begin from some initial guess for x(^), which is then refined via 
Hamilton's equations with the requirement that the 9i increase linearly with 
time. The two approaches are often complementary, as discussed below. 

2. Regular Motion 

In certain special potentials, every orbit is regular; examples are the Kepler 
and Stackel potentials. Motion in such potentials can be expressed most sim- 
ply by finding a canonical transformation to coordinates (p, q) for which the 
Hamiltonian is independent of q, H = H{p); among all such coordinates, one 
particularly simple choice is the action-angle variables {Ji,9i), in terms of which 
the equations of motion are 

Ji = constant, 

r)ff 

Oi = n,t + el ^^ = jj.^ i = i,...,iv (1) 

(Landau & Lifshitz 1976; Goldstein 1980). The trajectory x(J,^) is periodic in 
each of the angle variables 9i, which may be restricted to the range < ^j < 27r. 
The Ji define the cross-sectional areas of the torus while the 9i define position on 
the torus (Figure 1). These tori are sometimes called "invariant" since a phase 
point that lies on a torus at any time will remain on it forever. 

Most potentials are not integrable, but regular orbits may still exist; in- 
deed these are the orbits for which torus construction machinery is designed. 
One expects that for a regular orbit in a non-integrable potential, a canonical 
transformation (x, v) -^ (J,^) can be found such that 

J, = 0, e, = Vti, i = l,...,N. (2) 

However there is no guarantee that the full Hamiltonian will be expressible 
as a continuous function of the Ji. In general, the map (x,v) — > (J,^) will 
be different for each orbit and will not exist for those trajectories that do not 
respect N isolating integrals (although approximate maps, valid for some limited 
span of time, may be derived for weakly chaotic trajectories). 

The uniform translation of a regular orbit on its torus implies that the 
motion in any canonical coordinates (x, v) is quasi-periodic: 

x(t) = ^XA:(J)exp[i(/feJ7i + mfcri2 + nfef]3)t] , 

k 

w{t) = ^Yk{3)(iyiv[i{h^i + T^k^2 + nkVt^)t\, (3) 

k 

with {Ik^'^^ki'i^k) integers. The Fourier transform of x(t) or v(t) will therefore 
consist of a set of spikes at discrete frequencies uj^ = Ik^i + 'n^k^2 + ^fc^s 
that are linear combinations of the N fundamental frequencies J7j, with spectral 
amplitudes Xfc(J) and Vfc(J). 




Figure 1. Invariant torus defining the motion of a regular orbit in 
a two-dimensional potential. The torus is determined by the values 
of the actions Ji and J2; the position of the trajectory on the torus 
is defined by the angles 9i and 62, which increase linearly with time, 
6i = 0,it + 6^ . 



3. Trajectory-Following Approaches 

The most straightforward, and probably the most robust, approach to torus 
construction is via Fourier analysis of the numerically-integrated trajectories 
(Percival 1974; Boozer 1982; Binney & Spergel 1982, 1984; Kuo-Petravic et al. 
1983; Eaker et al. 1984; Martens & Ezra 1985). The Fourier decomposition 
of a quasiperiodic orbit (Equation 3) yields a discrete frequency spectrum. The 
precise form of this spectrum depends on the coordinates in which the orbit is in- 
tegrated, but certain of its properties are invariant, including the N fundamental 
frequencies fli from which every line is made up, uJk = /fcili+mfcri2 + f^fcf^3. Typ- 
ically the strongest line in a spectrum lies at one of the fundamental frequencies; 
once the ilj have been identified, the integer vectors (/fc, wifci'^fc) corresponding 
to every line cok are uniquely defined, to within computational uncertainties. 
Approximations to the actions may then be computed using Percival's (1974) 
formulae; e.g. the action associated with ^i in a 3 DOF system is 



Ji = ^ /fc (/fc$7i + mfcil2 + nk^s) |Xfc| 
k 



(4) 



and similarly for J2 and J3, upon replacing the first factor in the summation by 
nik and Uk respectively. Finally, the maps (9 -^ x) are obtained by making the 
substitution Ojt -^ 9i in the spectrum, e.g. 

x{t) = ^ Afc(J)exp[i(/fer2i + mfcil2 + "-fc^3)i] 

k 

= ^ Afc ( J) exp [i {lk9i + mk92 + ^1^6*3)] 
k 
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= x{ei,92,e3). (5) 

Trajectory following algorithms are easily automated; for instance, integer pro- 
gramming may be used to recover the vectors {Ikif^k^nk) (Valluri & Merritt 
1998). 

Binney & Spergel (1982) pioneered the use of trajectory-following algo- 
rithms for galactic potentials. They integrated orbits for a time T and computed 
discrete Fourier transforms, yielding spectra in which each frequency spike was 
represented by a peak with finite width ~ tt/T centered on uj^. They then fitted 
these peaks to the expected functional form X^ sin[(a; — ujj^)T]/{lo — ujk) using a 
least-squares algorithm. They were able to recover the fundamental frequencies 
in a 2 DOF potential with an accuracy of ~ 0.1% after ~ 25 orbital periods. Bin- 
ney & Spergel (1984) used equation (4) to construct the "action map" for orbits 
in a principal plane of the triaxial logarithmic potential. Carpintero & Aguilar 
(1998) and Copin, Zhao & de Zeeuw (this volume) applied similar algorithms to 
motion in 2- and 3 DOF potentials. 

The accuracy of Fourier transform methods can be greatly improved by 
multiplying the time series with a windowing function before transforming. The 
result is a reduction in the amplitude of the side lobes of each frequency peak 
at the expense of a broadening of the peaks; the amplitude measurements are 
then effectively decoupled from any errors in the determination of the frequen- 
cies. Laskar (1988, 1990) developed this idea into a set of tools, the "numerical 
analysis of fundamental frequencies" (NAFF), which he applied to the analy- 
sis of weakly chaotic motion in the solar system. Laskar's algorithm recovers 
the fundamental frequencies with an error that falls off as T~^ (Laskar 1996), 
compared with ~ T^^ in algorithms like Binney & Spergel's (1982). Even for 
modest integration times of ~ 10^ orbital periods, the NAFF algorithm is able 
to recover fundamental frequencies with accuracies of ~ 10~^ or better in many 
potentials. The result is a very precise representation of the torus (Figure 2). 

One drawback of trajectory- following algorithms is the need to extract a 
large number of terms in the frequency spectrum in cases where the time de- 
pendence of the integration variables is very different from that of the angles. 
This problem may be dealt with by expressing the numerically-integrated orbit 
in terms of a set of coordinates that are closer to the angle variables before 
computing the Fourier transform; for instance, tube orbits are most efficiently 
expressed in the canonically-conjugate Poincare variables (related to cylindrical 
coordinates, e.g. Papaphilippou & Laskar 1996). 

Trajectory-following algorithms also suffer from the fundamental limitation 
that they must follow the trajectory sufficiently far to see the longest periodicities 
of the orbit. In other words, the trajectory must adequately sample the surface 
of its invariant torus. Near a resonant torus, i.e. a torus for which the Oj satisfy 
a relation 

N 
1=1 

with the ai integers, trajectories fill their tori very slowly, necessitating long 
integration intervals. However even for near-resonant orbits, one can still effi- 
ciently recover the terms in the spectrum associated with the "faster" angles, as 
well as a reasonable approximation to the "slow" frequency f^a associated with 
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Figure 2. Construction of a 2 DOF, box-orbit torus in a Stackel 

potential using the NAFF trajectory-following algorithm, (a) The or- 
bit and its actions, computed using Equation (4) with kmax terms. 
Dashed lines show the exact Jj. (b) The map y{9i,62)', dashed con- 
tours correspond to negative values of y. A{kmax) is the RMS error 
in the reconstructed map, calculated using an equation similar to (5); 
A ~ /c-2 
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Figure 3. Resonant tori, (a) A two-dimensional torus, shown here as 
a square with identified edges. The plotted trajectory satisfies a 2 : 1 
resonance between the fundamental frequencies, 0,i — 2i72 = (e.g. a 
"banana"), (b) A three-dimensional torus, shown here as a cube with 
identified sides. The shaded region is covered densely by a resonant 
trajectory for which 2r2i + ^2 — 2^^ = 0. This trajectory is not closed, 
but it is restricted by the resonance condition to a two-dimensional 
subset of the torus. The orbit in configuration space is thin (Figure 4). 



libration around the resonant orbit (Merritt & Valluri 1999), and for many ap- 
plications these are sufficient. If the precise dependence of the map on ^3 is also 
needed, one possible approach (Papaphilippou & Laskar 1998) is to integrate 
the equations of motion using as a time step the period associated with one of 
the fast angles, thus eliminating it from the spectrum. 

Since Fourier techniques focus on the frequency domain, they are particu- 
larly well suited to identifying regions of phase space occupied by resonances. 
They are also ideal for studying the effect of resonances on the structure of phase 
space, even in cases where the full tori are difficult to reconstruct. Resonant tori 
are places where perturbation expansions of integrable systems break down, due 
to the "problem of small denominators" . In perturbed (non-integrable) poten- 
tials, one expects stable resonant tori to generate regions of regular motion and 
unstable resonant tori to give rise to chaotic regions. Algorithms like NAFF 
allow one to construct a "frequency map" of the phase space: a plot of the 
ratios of the fundamental frequencies (0,1/^^3,0,2/^3) for a large a set of or- 
bits selected from a uniform grid in initial condition space. Resonances appear 
on the frequency map as lines, either densely filled lines in the case of stable 
resonances, or gaps in the case of unstable resonances; the frequency map is 
effectively a representation of the Arnold web (Laskar 1993). Papaphilippou 
& Laskar (1996, 1998), Wachlin & Ferraz-Mello (1998) and Valluri & Merritt 
(1998) used frequency maps to study the effect of resonances on the structure 
of phase space in triaxial potentials. 




Figure 4. The surface in configuration space filled by a resonant, or 
"thin," box orbit in a triaxial potential (Merritt &: Valluri 1999). The 
order of the resonance is (2, 1, —2), as in Figure 3b. The surface was 
plotted by representing the spatial coordinates {x,y, z) parametrically 
in terms of the two angles that define position on the resonant three- 
torus. 



Precisely resonant orbits can be reconstructed using trajectory- following 
algorithms in a particularly straightforward way. A resonance has the effect of 
restricting an orbit to a subset of its torus, reducing the number of independent 
angle variables by one; thus a 2 DOF trajectory is reduced to a closed curve and 
a 3 DOF trajectory becomes a thin sheet (Born 1960; Goldstein 1980; Figure 3). 
The two frequencies defining motion on a resonant three-torus may be taken to 
be 



a 



in terms of which 
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This definition is not unique since the orbit is not closed. The motion in Carte- 
sian coordinates then becomes 

x(t) = ^ Xfc exp f (/fefJi + mfcrJ2 + n^ris) t 
k 



= ^Xfcexpi {-lkm3+nkmi)Qy + {-lkm2 + mkmi)Qy t 
k 

= J2 ^k exp i {ik'nil^ + nik'n^o^ ) t 

k 

= ^Xfcexpi(/fe'e«+mfc'e(2)^ 
k 

The result is a set of parametric expressions for the Cartesian coordinates 
{x,y,z) in terms of the angles {O^^'^O^"^') that define position on the two-torus 
(Figure 4). 

4. Iterative Approaches 

Iterative approaches to torus construction consist of finding successively better 
approximations to the map 6 ^ x given some initial guess x(0); canonical 
perturbation theory is a special case, and in fact iterative schemes often reduce 
to perturbative methods in appropriate limits. Iterative algorithms were first 
developed in the context of semi-classical quantization for computing energy 
levels of bound molecular systems, and they are still best suited to assigning 
energies to actions, H(J). Most of the other quantities of interest to galactic 
dynamicists - e.g. the fundamental frequencies Oj - are not recovered with high 
accuracy by these algorithms. Iterative schemes also tend to be numerically 
unstable unless the initial guess is close to the true solution. On the other hand, 
iterative algorithms can be more efficient than trajectory- following methods for 
orbits that are near resonance. 

Ratcliff, Chang & Schwarzschild (1984) pioneered iterative schemes in galac- 
tic dynamics. They noted that the equations of motion of a 2 DOF regular orbit, 

a; = -^-, 2/ = -^-, (10) 

ox ay 

can be written in the form 

If one specifies Qi and 0,2 and treats d^/dx and d^/dy as functions of the 9i, 
equations (11) can be viewed as nonlinear differential equations for x{9i, 62) and 
y(^i)^2)- Ratcliff et al. expressed the coordinates as Fourier series in the angle 
variables, 

x(0)=5]X„e-^. (12) 

n 

Substituting (12) into (11) gives 

^(n-0)2x„e^"-^ = V$ (13) 



where the right hand side is again understood to be a function of the angles. 
Ratchff et al. truncated the Fourier series after a finite number of terms and 
required equations (13) to be satisfied on a grid of points around the torus. 
They then solved for the X„ by iterating from an initial guess. Convergence 
was found to be possible if the initial guess was close to the exact solution. A 
similar algorithm was developed for recovering tori in the case that the actions, 
rather than the frequencies, are specified a priori. Guerra &; Ratcliff (1990) 
applied these algorithms to motion in the plane of rotation of a nonaxisymmetric 
potential. 

Another iterative approach to torus construction was developed by Chap- 
man, Garrett & Miller (1976) in the context of semiclassical quantum theory. 
One begins by dividing the Hamiltonian H into separable and non-separable 
parts Hq and Hi, then seeks a generating function S that maps the known tori 
of Hq into tori of H. For a generating function of the i<2-type (Goldstein 1980), 
one has 

J(^,J') = §, 0'{e,3') = §j (14) 

where (J, 9) and (J', 9') are the action-angle variables of Hq and H respectively. 
The generator S is determined, for a specified J', by substituting the first of 
equations (14) into the Hamiltonian and requiring the result to be independent 
of 9. One then arrives at H{J'). Chapman et al. showed that a sufficiently 
general form for S is 

S{9,J') = 9-3'-iY,Sn{JY'^''^, (15) 

where the first term is the identity transformation, and they evaluated a number 
of iterative schemes for finding the Sn- One such scheme was found to recover 
the results of first-order perturbation theory after a single iteration. McGill & 
Binney (1990) refined the Chapman et al. algorithm and applied it to 2 DOF 
motion in the axisymmetric logarithmic potential. 

The generating function approach is not naturally suited to deriving the 
other quantities of interest to galactic dynamicists. For instance, equation (14) 
gives 9'(9) as a derivative of S, but since S must be computed separately for 
every J' its derivative is likely to be ill-conditioned. Binney & Kumar (1993) 
and Kaasalainen & Binney (1994a) discussed two schemes for finding 9' {9); the 
first requires the solution of a formally infinite set of equations, while the lat- 
ter requires multiple integrations of the equations of motion for each torus - 
effectively a trajectory- following scheme. 

Kaasalainen & Binney (1994a) noted that the success of the generating 
function method depends strongly on the choice of Hq. For box orbits, which 
are most naturally described as coupled rectilinear oscillators, they found that 
a harmonic-oscillator Hq gave poor results unless an additional point transfor- 
mation was used to deform the rectangular orbits of Hq into narrow-waisted 
boxes like those in typical galactic potentials. Kaasalainen (1995a) considered 
orbits belonging to higher-order resonant families and found that it was generally 
necessary to define a new coordinate transformation for each family. 

Warnock (1991) presented a hybrid scheme in which the generating func- 
tion S was derived by numerically integrating an orbit from appropriate initial 
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conditions, transforming the coordinates to (J, 6) of Hq and interpolating J on 
a regular grid in 9. The values of the Sn then follow from the first equation of 
(14) after a discrete Fourier transform. Kaasalainen & Binney (1994b) found 
that Warnock's scheme could be used to substantially refine the solutions found 
via their iterative algorithm. Another hybrid scheme was discussed by Reiman 
& Pomphrey (1991). 

Having computed the energy on a grid of J' values, one can interpolate 
to obtain the full Hamiltonian H{J'). If the system is not in fact completely 
integrable, this H may be rigorously interpreted as smooth approximation to 
the true H (Warnock & Ruth 1991, 1992) and can be taken as the starting 
point for secular perturbation theory. Kaasalainen (1994) developed this idea 
and showed how to recover accurate surfaces of section in the neighborhood of 
low-order resonances in the planar logarithmic potential. 

Percival (1977) described a variational principle for constructing tori. His 
technique has apparently not yet been implemented in the context of galactic 
dynamics. 

5. Chaotic Motion 

Torus-construction machinery may be applied to orbits that are approximately, 
but not precisely, regular (Laskar 1993). The frequency spectrum of a weakly 
chaotic orbit will typically be close to that of a regular orbit, with most of the 
lines well approximated as linear combinations of three "fundamental frequen- 
cies" Qi. However these frequencies will change with time as the orbit diffuses 
from one "torus" to another. The diffusion rate can be measured via quantities 
like I ill — fi'il, the change in a "fundamental frequency" over two consecutive 
integration intervals. Papaphilippou &: Laskar (1996, 1998), Valluri & Merritt 
(1998) and Wachlin & Ferraz-Mello (1998) used this technique to study chaos 
and diffusion in triaxial galactic potentials. Kaasalainen (1995b) showed that 
approximate tori could be constructed even in chaotic phase space via the hybrid 
scheme of Warnock (1991). While such tori clearly do not describe the motion 
of chaotic orbits over long times, they are useful for understanding the onset of 
chaos and its relationship to resonances, as well as for studying evolution of the 
phase-space distribution function in action space via the Fokker-Plank equation 
(Lichtenberg & Leiberman 1992). 

6. Summary 

Trajectory-following schemes for torus construction are robust and easily auto- 
mated. They can recover the fundamental frequencies with great precision and 
are well suited to studies of weak chaos and for mapping resonances. However 
they are inefficient for constructing the full torus of an orbit that lies close to, 
but slightly off of, a resonance. Iterative techniques are efficient for assigning 
energies to actions but less suited to recovering the other quantities of interest 
to galactic dynamicists, such as the fundamental frequencies. However they can 
be more efficient than trajectory-following algorithms for constructing nearly- 
resonant tori. Hybrid schemes that combine features of both approaches show 
considerable promise. 
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